#ifndef _SPLINE_H_
#define _SPLINE_H_

#include "Polynomial.h"
#include "Numtool.h"
#include "InterpConditions.h"
#include "Bspline.h"
#include "KnotsHelper.h"
#include <fstream>
#include <string>
#include <cblas.h>
#include <lapacke.h>

enum SplineType{
            ppForm,
            cardinalB,
            BForm         
};

enum BCType{
              complete,
              notAknot,
              periodic,
              second,
              mixed,
              meaningless
};

template<int Dim, int Order, SplineType t>
class Spline;

template<int Order,int NUM,SplineType t>
Spline<1,Order,t> interpolate(const InterpConditions<double,NUM>& IC,
                              BCType type = meaningless,
                              std::pair<int,int> infopair = std::make_pair(0,0));

template<int Dim, int Order>
class Spline<Dim,Order,ppForm>{
public:
  using PolyDim = Polynomial<Order,Vec<double,Dim> >;
protected:
  std::vector<double> knots_;
  std::vector<PolyDim> polys_;
public:
  Spline(const double startknot = 0);
  Spline(const PolyDim& p, const double leng, const double startknot =
         0);
  ~Spline() = default;
  const std::vector<double>& getknots() const;
  void setknots(const std::vector<double>& knt);
  const std::vector<Polynomial<Order,Vec<double,Dim> > >& getpolys() const;
  void setpolys(const std::vector<PolyDim>& p);
  void addknot(const PolyDim& p, const double leng);
  void addcurve(const Spline<Dim,Order,ppForm>& ppsp);
  int locatePiece(double t) const;
  Vec<double,Dim> operator()(double t) const;
  template<int Ord>
  friend Spline<1,Ord,ppForm> operator+(const Spline<1,Ord,ppForm>& sp1,
                                        const Spline<1,Ord,ppForm>& sp2);
  template<int Ord>
  friend Spline<1,Ord,ppForm> operator-(const Spline<1,Ord,ppForm>& sp1,
                                        const Spline<1,Ord,ppForm>& sp2);
  void show(std::ostream& os = std::cout) const;
  template<int Ord>
  friend void fnplt(const Spline<1,Ord,ppForm>& sp, const std::string
                    &_file, const int num, const bool overwrite);
  template<int Ord>
  friend void fnplt(const Spline<2,Ord,ppForm>& sp, const std::string
                    &_file, const int num, const bool overwrite, const bool isperiod);
  template<int D,int Ord,int InputD>
  friend Spline<D,Ord,ppForm> fitCurve(const
                                       std::vector<Vec<double,InputD>
                                       >& nknots,
                                       BCType type,
                                       const Vec<double,D>& start,
                                       const Vec<double,D>& end,
                                       std::pair<int,int> infopair);
  template<int Ord,int num,SplineType T>
  friend Spline<1,Ord,T> interpolate(const
                                     InterpConditions<double,num>& IC,
                                     BCType type, std::pair<int,int> infopair);
};

template<int Dim, int Order>
inline Spline<Dim,Order,ppForm>::Spline(const double startknot){
  knots_.push_back(startknot);
}

template<int Dim, int Order>
inline Spline<Dim,Order,ppForm>::Spline(const PolyDim& p, const double leng, const double startknot){
  knots_.push_back(startknot);
  knots_.push_back(startknot+leng);
  polys_.push_back(p);
}

template<int Dim, int Order>
inline const std::vector<double>& Spline<Dim,Order,ppForm>::getknots() const{
  return knots_;
}

template<int Dim, int Order>
inline void Spline<Dim,Order,ppForm>::setknots(const std::vector<double>& knt){
  knots_ = knt;
}


template<int Dim, int Order>
inline const std::vector<Polynomial<Order,Vec<double,Dim> > >& Spline<Dim,Order,ppForm>::getpolys() const{
  return polys_;
}

template<int Dim, int Order>
inline void Spline<Dim,Order,ppForm>::setpolys(const std::vector<PolyDim>& p){
  polys_ = p;
}

template<int Dim, int Order>
inline void Spline<Dim,Order,ppForm>::addknot(const PolyDim& p, const double leng){
  assert(!knots_.empty());
  knots_.push_back(knots_.back()+leng);
  polys_.push_back(p);
}

template<int Dim, int Order>
void Spline<Dim,Order,ppForm>::addcurve(const Spline<Dim,Order,ppForm>& ppsp){
  assert(!knots_.empty());
  const std::vector<PolyDim>& sppolys = ppsp.getpolys();
  if (sppolys.empty())
    return;
  const std::vector<double>& spknots = ppsp.getknots();
  const int Size = sppolys.size();
  const double d = knots_.back() - spknots.front();
  for (int i = 0; i < Size ; i++){
    knots_.push_back(d+spknots[i+1]);
    polys_.push_back(sppolys[i]);
  }    
}

template<int Dim, int Order>
int Spline<Dim,Order,ppForm>::locatePiece(double t) const{
  constexpr double tol = 1e-10;
  assert(t > knots_.front() - tol && t < knots_.back() + tol);
  const int num = knots_.size();
  int left = 0, right = num - 2, mid;
  while(1){
    if (left == right)
      return left;
    mid = (left+right)/2;
    if (t < knots_[mid+1])
      right = mid;
    else
      left = mid+1;
  }
}

template<int Dim, int Order>
inline Vec<double,Dim> Spline<Dim,Order,ppForm>::operator()(double t) const{
  int index = locatePiece(t);
  return (polys_[index])(t-knots_[index]);
}



template<int Order>
Spline<1,Order,ppForm> operator+(const Spline<1,Order,ppForm>& sp1,
                                 const Spline<1,Order,ppForm>& sp2){
  const std::vector<double>& knots1 = sp1.getknots();
  const int N1 = knots1.size();
  const std::vector<double>& knots2 = sp2.getknots();
  const int N2 = knots2.size();
  std::vector<double> knots3(N1+N2);
  std::merge(knots1.cbegin(),knots1.cend(),knots2.cbegin(),knots2.cend(),knots3.begin());
  auto pos = std::unique(knots3.begin(),knots3.end());
  knots3.erase(pos,knots3.end());
  const int num = knots3.size();
  const std::vector<Polynomial<Order,Vec<double,1> > >& polys1 =
  sp1.getpolys();
  const std::vector<Polynomial<Order,Vec<double,1> > >& polys2 =
  sp2.getpolys();
  Spline<1,Order,ppForm> res;
  res.setknots(knots3);
  std::vector<Polynomial<Order,Vec<double,1> > > respoly(num-1);
  for (int i = 0; i < num - 1 ; i++){
    double tmp = knots3[i];
    int d1 = sp1.locatePiece(tmp);
    // if (equal(knots3[i],knots1[d1+1]))
    //   d1++;
    int d2 = sp2.locatePiece(tmp);
    // if (equal(knots3[i],knots2[d2+1]))
    //   d2++;
    respoly[i] = (polys1[d1].translate(tmp - knots1[d1]) +
    polys2[d2].translate(tmp - knots2[d2]));
  }
  res.setpolys(respoly);
  return res;
}

template<int Order>
Spline<1,Order,ppForm> operator-(const Spline<1,Order,ppForm>& sp1,
                                 const Spline<1,Order,ppForm>& sp2){
  const std::vector<double>& knots1 = sp1.getknots();
  const int N1 = knots1.size();
  const std::vector<double>& knots2 = sp2.getknots();
  const int N2 = knots2.size();
  std::vector<double> knots3(N1+N2);
  std::merge(knots1.cbegin(),knots1.cend(),knots2.cbegin(),knots2.cend(),knots3.begin());
  auto pos = std::unique(knots3.begin(),knots3.end());
  knots3.erase(pos,knots3.end());
  const int num = knots3.size();
  const std::vector<Polynomial<Order,Vec<double,1> > >& polys1 =
  sp1.getpolys();
  const std::vector<Polynomial<Order,Vec<double,1> > >& polys2 =
  sp2.getpolys();
  Spline<1,Order,ppForm> res;
  res.setknots(knots3);
  std::vector<Polynomial<Order,Vec<double,1> > > respoly(num-1);
  for (int i = 0; i < num - 1 ; i++){
    double tmp = knots3[i];
    int d1 = sp1.locatePiece(tmp);
    // if (equal(tmp,knots1[d1+1]))
    //   d1++;
    int d2 = sp2.locatePiece(tmp);
    // if (equal(tmp,knots2[d2+1]))
    //   d2++;
    respoly[i] = (polys1[d1].translate(tmp - knots1[d1]) -
    polys2[d2].translate(tmp - knots2[d2]));
  }
  res.setpolys(respoly);
  return res;
}


template<int Dim, int Order>
inline void Spline<Dim,Order,ppForm>::show(std::ostream& os) const{
  for (int i = 0 ; i < (int)polys_.size() ; i++)
    os << "[" << knots_[i] << "," << knots_[i+1] << "]: " << polys_[i] << std::endl;
}

template<int Order>
void fnplt(const Spline<1,Order,ppForm>& sp, const std::string
           &_file, const int num = 200, const bool overwrite = 0){
  const std::vector<Polynomial<Order,Vec<double,1> > >& sppolys =
  sp.getpolys();
  assert(!sppolys.empty());
  const std::vector<double>& spknots = sp.getknots();
  std::ofstream os;
  if (overwrite)
    os.open(_file,ios::app);
  else
    os.open(_file);
  double h = (spknots.back() - spknots.front())/num;
  std::vector<double> y(num+1);
  int j = 0;
  double t = spknots.front();
  os << "t = [";
  for (int i = 0; i < num ; i++){
    if (t >= spknots[j+1])
      j++;
    Vec<double,1> tmpvec = (sppolys[j])(t-spknots[j]);
    y[i] = tmpvec[0];
    os << t << ",";
    t += h;
  }
  os << t;
  Vec<double,1> tmpvec = (sppolys[j])(t-spknots[j]);
  y[num] = tmpvec[0];
  os << "];\n" << "y = [";
  for (int i = 0; i <= num ; i++)
    os << y[i] << ",";
  os << "];\n" << "plot(t,y);";
  os.close();
}


template<int Order>
void fnplt(const Spline<2,Order,ppForm>& sp, const std::string &_file,
           const int num = 200,  const bool
           overwrite = 0, const bool isperiod = 0){
  const std::vector<Polynomial<Order,Vec<double,2> > >& sppolys = sp.getpolys();
  assert(!sppolys.empty());
  const std::vector<double>& spknots = sp.getknots();
  std::ofstream os;
  if (overwrite)
    os.open(_file,ios::app);
  else
    os.open(_file);
  double h = (spknots.back() - spknots.front())/num;
  std::vector<double> x(num+1),y(num+1);
  int j = 0;
  double t = spknots.front();
  os << "t = [";
  for (int i = 0; i < num ; i++){
    if (t >= spknots[j+1])
      j++;
    Vec<double,2> tmpvec = (sppolys[j])(t-spknots[j]);
    x[i] = tmpvec[0];
    y[i] = tmpvec[1];
    os << t << ",";
    t += h;
  }
  os << t;
  Vec<double,2> tmpvec = (sppolys[j])(t-spknots[j]);
  x[num] = tmpvec[0];
  y[num] = tmpvec[1];
  os << "];\n" << "x = [";
  for (int i = 0; i < num ; i++)
    os << x[i] << ",";
  if (isperiod)
    os << x[0];
  else
    os << x[num];
  os << "];\n" << "y = [";
  for (int i = 0; i < num ; i++)
    os << y[i] << ",";
  if (isperiod)
    os << y[0];
  else
    os << y[num];
  os << "];\n" << "plot(x,y);";
  os.close();
}


template<int Dim,int Order,int InputD = Max(2,Dim)>
Spline<Dim,Order,ppForm> fitCurve(const std::vector<Vec<double,InputD>
                                   >& nknots,
                                   BCType type = meaningless,
                                   const Vec<double,Dim>& start = Vec<double,Dim>(),
                                   const Vec<double,Dim>& end = Vec<double,Dim>(),
                                   std::pair<int,int> infopair= std::pair<int,int>());



template<> Spline<2,2,ppForm> fitCurve(const
                                   std::vector<Vec<double,2>>& nknots,
                                   BCType type,
                                   const Vec<double,2>& start,
                                   const Vec<double,2>& end,
                                   std::pair<int,int> infopair){
  assert(type == meaningless);
  const int Size = nknots.size();
  assert(Size >= 2);
  Spline<2,2,ppForm> res;
  for (int i = 0 ; i < Size - 1 ; i++){
    double l = norm(nknots[i+1] - nknots[i],2);
    Polynomial<2,Vec<double,2> > p;
    p[0] = (nknots[i+1] - nknots[i])/l;
    p[1] = nknots[i];
    res.addknot(p,l);
  }
  return res;
}


template<> Spline<2,4,ppForm> fitCurve(const
                                   std::vector<Vec<double,2>>& nknots,
                                   BCType type,
                                   const Vec<double,2>& start,
                                   const Vec<double,2>& end,
                                   std::pair<int,int> infopair){
  if (type == meaningless) //default case
    type = notAknot;
  const int N = nknots.size() - 1;
  assert(N >= 1);
  std::vector<double> t(N+1);
  t[0] = 0.0;
  for (int i = 1; i <= N ; i++)
    t[i] = t[i-1] + norm(nknots[i]-nknots[i-1],2);
  const auto& v = nknots;
  Spline<2,4,ppForm> res;
  //periodic case:
  if (type == periodic){
    //coefficient matrix
    assert(equal<2>(nknots[0],nknots[N]));
    double* M = new double[N*N];
    for (int i = 0 ; i < N*N ; i++)
      M[i] = 0.0;
    M[0] = 2.0;
    M[N] = (t[N] - t[N-1])/(t[1] + t[N] - t[N-1]);
    M[N*(N-1)] = t[1]/(t[1] + t[N] - t[N-1]);
    for (int i = 1 ; i < N - 1 ; i++){
      M[i*N+i] = 2.0;
      M[(i-1)*N+i] = (t[i+1] - t[i])/(t[i+1] - t[i-1]);
      M[(i+1)*N+i] = (t[i] - t[i-1])/(t[i+1] - t[i-1]);
    }
    M[N*N-1] = 2.0;
    M[N*(N-1)-1] = (t[N] - t[N-1])/(t[N] - t[N-2]);
    M[N-1] = (t[N-1] - t[N-2])/(t[N] - t[N-2]);
    //prepare RHS
    double* RHS = new double[N*2];
    RHS[0] = 3*((t[N] - t[N-1])/(t[1] + t[N] - t[N-1]))*((v[1][0] -
                                   v[0][0])/t[1])+ 3*(t[1]/(t[1] +
                                   t[N] - t[N-1]))*((v[0][0] -
                                   v[N-1][0])/(t[N] - t[N-1]));
    RHS[N] = 3*((t[N] - t[N-1])/(t[1] + t[N] - t[N-1]))*((v[1][1] -
                                   v[0][1])/t[1])+ 3*(t[1]/(t[1] +
                                   t[N] - t[N-1]))*((v[0][1] -
                                   v[N-1][1])/(t[N] - t[N-1]));
    for (int i = 1 ; i < N ; i++){
      RHS[i] = 3*((t[i] - t[i-1])/(t[i+1] -
                                   t[i-1]))*((v[i+1][0]-v[i][0])/(t[i+1]-t[i]))+3*((t[i+1]
                                   - t[i])/(t[i+1] -
                                   t[i-1]))*((v[i][0]-v[i-1][0])/(t[i]-t[i-1]));
      RHS[N+i] = 3*((t[i] - t[i-1])/(t[i+1] -
                                   t[i-1]))*((v[i+1][1]-v[i][1])/(t[i+1]-t[i]))+3*((t[i+1]
                                   - t[i])/(t[i+1] -
                                   t[i-1]))*((v[i][1]-v[i-1][1])/(t[i]-t[i-1]));
    }
    //solve the solution
    int* ipiv = new int[N];
    auto info = LAPACKE_dgesv(LAPACK_COL_MAJOR,N,2,M,N,ipiv,RHS,N);
    if (info != 0)
      throw std::runtime_error("fitCurve<2,4,ppForm>() - DGESV");
    //assemble the spline
    res.knots_ = t;
    res.polys_.resize(N);
    for (int i = 0 ; i < N - 1 ; i++){
      Vec<double,2> K({(v[i+1][0]-v[i][0])/(t[i+1]-t[i]),(v[i+1][1]-v[i][1])/(t[i+1]-t[i])});
      Vec<double,2> m1({RHS[i],RHS[N+i]});
      Vec<double,2> m2({RHS[i+1],RHS[N+i+1]});
      res.polys_[i][0] = (m1 + m2 - (K * 2))/((t[i+1] - t[i])*(t[i+1] -
                                   t[i]));
      res.polys_[i][1] = (K * 3 - (m1 * 2) - m2)/(t[i+1] - t[i]);
      res.polys_[i][2] = m1;
      res.polys_[i][3] = Vec<double,2> {v[i][0],v[i][1]};
    }
    Vec<double,2> K({(v[N][0]-v[N-1][0])/(t[N]-t[N-1]),(v[N][1]-v[N-1][1])/(t[N]-t[N-1])});
    Vec<double,2> m1({RHS[N-1],RHS[2*N-1]});
    Vec<double,2> m2({RHS[0],RHS[N]});
    res.polys_[N-1][0] = (m1 + m2 - (K * 2))/((t[N] - t[N-1])*(t[N] -
                                   t[N-1]));
    res.polys_[N-1][1] = (K * 3 - (m1 * 2) - m2)/(t[N] - t[N-1]);
    res.polys_[N-1][2] = m1;
    res.polys_[N-1][3] = Vec<double,2> {v[N-1][0],v[N-1][1]};
    //release memory
    delete [] M;
    delete [] RHS;
    delete [] ipiv;
  }
  //complete,second,notAknot and variational case:
  if (type == complete || type == second || type == notAknot || type
                                   == mixed){
    //determine conditions
    if (type == complete)
      infopair = std::make_pair(1,1);
    else if (type == second)
      infopair = std::make_pair(2,2);
    else if (type == notAknot)
      infopair = std::make_pair(0,0);
    //coefficient matrix
    double* M = new double[(N+1)*(N+1)];
    for (int i = 0 ; i < (N+1)*(N+1); i++)
      M[i] = 0.0;
    if (infopair.first == 1)
      M[0] = 1.0;
    else if (infopair.first == 2){
      M[0] = 2.0;
      M[N+1] = 1.0;
    }
    else if (infopair.first == 0){
      double tmp = (t[1] - t[0])*(t[1] - t[0])/((t[2] - t[1])*(t[2] - t[1]));
      M[0] = 1.0;
      M[N+1] = 1 - tmp;
      M[2*N+2] = -tmp;  
    }
    else{
      throw std::runtime_error("First parameter of infopair is out of \
                                range! It should be 0,1 or 2.");
    }
    for (int i = 1 ; i < N ; i++){
      M[i*(N+1)+i] = 2.0;
      M[(i-1)*(N+1)+i] = (t[i+1] - t[i])/(t[i+1] - t[i-1]);
      M[(i+1)*(N+1)+i] = (t[i] - t[i-1])/(t[i+1] - t[i-1]);
    }
    if (infopair.second == 1)
      M[(N+1)*(N+1)-1] = 1.0;
    else if (infopair.second == 2){
      M[(N+1)*(N+1)-1] = 2.0;
      M[N*(N+1)-1] = 1.0;
    }
    else if (infopair.second == 0){
      double tmp = (t[N] - t[N-1])*(t[N] - t[N-1])/((t[N-1] - t[N-2])*(t[N-1] - t[N-2]));
      M[(N+1)*(N+1)-1] = 1.0;
      M[N*(N+1)-1] = 1 - tmp;
      M[(N-1)*(N+1)-1] = -tmp;  
    }
    else{
      throw std::runtime_error("Second parameter of infopair is out of\
                                   range! It should be 0,1 or 2.");
    }
    //prepare RHS
    double* RHS = new double[(N+1)*2];
    if (infopair.first == 1){
      RHS[0] = start[0];
      RHS[N+1] = start[1];
    }
    else if (infopair.first == 2){
      RHS[0] = 3*((v[1][0]-v[0][0])/(t[1]-t[0]))-0.5*start[0]*(t[1]-t[0]);
      RHS[N+1] = 3*((v[1][1]-v[0][1])/(t[1]-t[0]))-0.5*start[1]*(t[1]-t[0]);
    }
    else if (infopair.first == 0){
      double tmp = (t[1] - t[0])*(t[1] - t[0])/((t[2] - t[1])*(t[2] -
                                   t[1]));
      RHS[0] = -2*tmp*((v[2][0]-v[1][0])/(t[2]-t[1]))+2*((v[1][0]-v[0][0])/(t[1]-t[0]));
      RHS[N+1] = -2*tmp*((v[2][1]-v[1][1])/(t[2]-t[1]))+2*((v[1][1]-v[0][1])/(t[1]-t[0]));
    }
    for (int i = 1 ; i < N; i++){
      RHS[i] = 3*((t[i] - t[i-1])/(t[i+1] -
                                   t[i-1]))*((v[i+1][0]-v[i][0])/(t[i+1]-t[i]))+3*((t[i+1]
                                   - t[i])/(t[i+1] -
                                   t[i-1]))*((v[i][0]-v[i-1][0])/(t[i]-t[i-1]));
      RHS[N+i+1] = 3*((t[i] - t[i-1])/(t[i+1] -
                                   t[i-1]))*((v[i+1][1]-v[i][1])/(t[i+1]-t[i]))+3*((t[i+1]
                                   - t[i])/(t[i+1] -
                                   t[i-1]))*((v[i][1]-v[i-1][1])/(t[i]-t[i-1]));
    }
    if (infopair.second == 1){
      RHS[N] = end[0];
      RHS[2*N+1] = end[1];
    }
    else if (infopair.second == 2){
      RHS[N] = 3*((v[N][0]-v[N-1][0])/(t[N]-t[N-1]))-0.5*end[0]*(t[N]-t[N-1]);
      RHS[2*N+1] = 3*((v[N][1]-v[N-1][1])/(t[N]-t[N-1]))-0.5*end[1]*(t[N]-t[N-1]);
    }
    else if (infopair.second == 0){
      double tmp = (t[N] - t[N-1])*(t[N] - t[N-1])/((t[N-1] -
                                   t[N-2])*(t[N-1] - t[N-2]));
      RHS[N] = -2*tmp*((v[N-1][0]-v[N-2][0])/(t[N-1]-t[N-2]))+2*((v[N][0]-v[N-1][0])/(t[N]-t[N-1]));
      RHS[2*N+1] = -2*tmp*((v[N-1][1]-v[N-2][1])/(t[N-1]-t[N-2]))+2*((v[N][1]-v[N-1][1])/(t[N]-t[N-1]));
    }
    //solve the solution
    int* ipiv = new int[N+1];
    auto info = LAPACKE_dgesv(LAPACK_COL_MAJOR,N+1,2,M,N+1,ipiv,RHS,N+1);
    if (info != 0)
      throw std::runtime_error("fitCurve<2,4,ppForm>() - DGESV");
    //assemble the spline
    res.knots_ = t;
    res.polys_.resize(N);
    for (int i = 0 ; i < N ; i++){
      Vec<double,2> K({(v[i+1][0]-v[i][0])/(t[i+1]-t[i]),(v[i+1][1]-v[i][1])/(t[i+1]-t[i])});
      Vec<double,2> m1({RHS[i],RHS[N+i+1]});
      Vec<double,2> m2({RHS[i+1],RHS[N+i+2]});
      res.polys_[i][0] = (m1 + m2 - (K * 2))/((t[i+1] - t[i])*(t[i+1] -
                                   t[i]));
      res.polys_[i][1] = (K * 3 - (m1 * 2) - m2)/(t[i+1] - t[i]);
      res.polys_[i][2] = m1;
      res.polys_[i][3] = Vec<double,2> {v[i][0],v[i][1]};
    }
    //release memory
    delete [] M;
    delete [] RHS;
    delete [] ipiv;
  }
  return res;
}

template<> Spline<1,2,ppForm> fitCurve(const
                                   std::vector<Vec<double,2>>& nknots,
                                   BCType type,
                                   const Vec<double,1>& start,
                                   const Vec<double,1>& end,
                                   std::pair<int,int> infopair){
  assert(type == meaningless);
  const int Size = nknots.size();
  assert(Size >= 2);
  Spline<1,2,ppForm> res(nknots[0][0]);
  for (int i = 0 ; i < Size - 1 ; i++){
    double l = nknots[i+1][0] - nknots[i][0];
    Polynomial<2,Vec<double,1> > p;
    p[0] = (nknots[i+1][1] - nknots[i][1])/l;
    p[1] = nknots[i][1];
    res.addknot(p,l);
  }
  return res;
}

template<> Spline<1,4,ppForm> fitCurve(const
                                   std::vector<Vec<double,2>>& nknots,
                                   BCType type,
                                   const Vec<double,1>& start,
                                   const Vec<double,1>& end,
                                   std::pair<int,int> infopair){
  if (type == meaningless) //default case
    type = notAknot;
  const int N = nknots.size() - 1;
  assert(N >= 1);
  std::vector<double> t(N+1);
  t[0] = nknots[0][0];
  for (int i = 1; i <= N ; i++)
    t[i] = t[i-1] + nknots[i][0] - nknots[i-1][0];
  const auto& v = nknots;
  Spline<1,4,ppForm> res;
  //periodic case:
  if (type == periodic){
    //coefficient matrix
    assert(equal<2>(nknots[0][1],nknots[N][1]));
    double* M = new double[N*N];
    for (int i = 0 ; i < N*N ; i++)
      M[i] = 0.0;
    M[0] = 2.0;
    M[N] = (t[N] - t[N-1])/(t[1] + t[N] - t[N-1]);
    M[N*(N-1)] = t[1]/(t[1] + t[N] - t[N-1]);
    for (int i = 1 ; i < N - 1 ; i++){
      M[i*N+i] = 2.0;
      M[(i-1)*N+i] = (t[i+1] - t[i])/(t[i+1] - t[i-1]);
      M[(i+1)*N+i] = (t[i] - t[i-1])/(t[i+1] - t[i-1]);
    }
    M[N*N-1] = 2.0;
    M[N*(N-1)-1] = (t[N] - t[N-1])/(t[N] - t[N-2]);
    M[N-1] = (t[N-1] - t[N-2])/(t[N] - t[N-2]);
    //prepare RHS
    double* RHS = new double[N];
    RHS[0] = 3*((t[N] - t[N-1])/(t[1] + t[N] - t[N-1]))*((v[1][1] -
                                   v[0][1])/t[1])+ 3*(t[1]/(t[1] +
                                   t[N] - t[N-1]))*((v[0][1] -
                                   v[N-1][1])/(t[N] - t[N-1]));
    for (int i = 1 ; i < N ; i++){
      RHS[i] = 3*((t[i] - t[i-1])/(t[i+1] -
                                   t[i-1]))*((v[i+1][1]-v[i][1])/(t[i+1]-t[i]))+3*((t[i+1]
                                   - t[i])/(t[i+1] -
                                   t[i-1]))*((v[i][1]-v[i-1][1])/(t[i]-t[i-1]));
    }
    //solve the solution
    int* ipiv = new int[N];
    auto info = LAPACKE_dgesv(LAPACK_COL_MAJOR,N,1,M,N,ipiv,RHS,N);
    if (info != 0)
      throw std::runtime_error("fitCurve<1,4,ppForm>() - DGESV");
    //assemble the spline
    res.knots_ = t;
    res.polys_.resize(N);
    for (int i = 0 ; i < N - 1 ; i++){
      Vec<double,1> K({(v[i+1][1]-v[i][1])/(t[i+1]-t[i])});
      Vec<double,1> m1({RHS[i]});
      Vec<double,1> m2({RHS[i+1]});
      res.polys_[i][0] = (m1 + m2 - (K * 2))/((t[i+1] - t[i])*(t[i+1] -
                                   t[i]));
      res.polys_[i][1] = (K * 3 - (m1 * 2) - m2)/(t[i+1] - t[i]);
      res.polys_[i][2] = m1;
      res.polys_[i][3] = Vec<double,1> {v[i][1]};
    }
    Vec<double,1> K({(v[N][1]-v[N-1][1])/(t[N]-t[N-1])});
    Vec<double,1> m1({RHS[N-1]});
    Vec<double,1> m2({RHS[0]});
    res.polys_[N-1][0] = (m1 + m2 - (K * 2))/((t[N] - t[N-1])*(t[N] -
                                   t[N-1]));
    res.polys_[N-1][1] = (K * 3 - (m1 * 2) - m2)/(t[N] - t[N-1]);
    res.polys_[N-1][2] = m1;
    res.polys_[N-1][3] = Vec<double,1> {v[N-1][1]};
    //release memory
    delete [] M;
    delete [] RHS;
    delete [] ipiv;
  }
  //complete,second,notAknot and variational case:
  if (type == complete || type == second || type == notAknot || type
                                   == mixed){
    //determine conditions
    if (type == complete)
      infopair = std::make_pair(1,1);
    else if (type == second)
      infopair = std::make_pair(2,2);
    else if (type == notAknot)
      infopair = std::make_pair(0,0);
    //coefficient matrix
    double* M = new double[(N+1)*(N+1)];
    for (int i = 0 ; i < (N+1)*(N+1); i++)
      M[i] = 0.0;
    if (infopair.first == 1)
      M[0] = 1.0;
    else if (infopair.first == 2){
      M[0] = 2.0;
      M[N+1] = 1.0;
    }
    else if (infopair.first == 0){
      double tmp = (t[1] - t[0])*(t[1] - t[0])/((t[2] - t[1])*(t[2] - t[1]));
      M[0] = 1.0;
      M[N+1] = 1 - tmp;
      M[2*N+2] = -tmp;  
    }
    else{
      throw std::runtime_error("First parameter of infopair is out of \
                                  range! It should be 0,1 or 2.");
    }
    for (int i = 1 ; i < N ; i++){
      M[i*(N+1)+i] = 2.0;
      M[(i-1)*(N+1)+i] = (t[i+1] - t[i])/(t[i+1] - t[i-1]);
      M[(i+1)*(N+1)+i] = (t[i] - t[i-1])/(t[i+1] - t[i-1]);
    }
    if (infopair.second == 1)
      M[(N+1)*(N+1)-1] = 1.0;
    else if (infopair.second == 2){
      M[(N+1)*(N+1)-1] = 2.0;
      M[N*(N+1)-1] = 1.0;
    }
    else if (infopair.second == 0){
      double tmp = (t[N] - t[N-1])*(t[N] - t[N-1])/((t[N-1] - t[N-2])*(t[N-1] - t[N-2]));
      M[(N+1)*(N+1)-1] = 1.0;
      M[N*(N+1)-1] = 1 - tmp;
      M[(N-1)*(N+1)-1] = -tmp;  
    }
    else{
      throw std::runtime_error("Second parameter of infopair is out of\
                                     range! It should be 0,1 or 2.");
    }
    //prepare RHS
    double* RHS = new double[N+1];
    if (infopair.first == 1){
      RHS[0] = start[0];
    }
    else if (infopair.first == 2){
      RHS[0] = 3*((v[1][1]-v[0][1])/(t[1]-t[0]))-0.5*start[0]*(t[1]-t[0]);
    }
    else if (infopair.first == 0){
      double tmp = (t[1] - t[0])*(t[1] - t[0])/((t[2] - t[1])*(t[2] -
                                   t[1]));
      RHS[0] = -2*tmp*((v[2][1]-v[1][1])/(t[2]-t[1]))+2*((v[1][1]-v[0][1])/(t[1]-t[0]));
    }
    for (int i = 1 ; i < N; i++){
      RHS[i] = 3*((t[i] - t[i-1])/(t[i+1] -
                                   t[i-1]))*((v[i+1][1]-v[i][1])/(t[i+1]-t[i]))+3*((t[i+1]
                                   - t[i])/(t[i+1] -
                                   t[i-1]))*((v[i][1]-v[i-1][1])/(t[i]-t[i-1]));
    }
    if (infopair.second == 1){
      RHS[N] = end[0];
    }
    else if (infopair.second == 2){
      RHS[N] = 3*((v[N][1]-v[N-1][1])/(t[N]-t[N-1]))-0.5*end[0]*(t[N]-t[N-1]);
    }
    else if (infopair.second == 0){
      double tmp = (t[N] - t[N-1])*(t[N] - t[N-1])/((t[N-1] -
                                   t[N-2])*(t[N-1] - t[N-2]));
      RHS[N] = -2*tmp*((v[N-1][1]-v[N-2][1])/(t[N-1]-t[N-2]))+2*((v[N][1]-v[N-1][1])/(t[N]-t[N-1]));
    }
    //solve the solution
    int* ipiv = new int[N+1];
    auto info = LAPACKE_dgesv(LAPACK_COL_MAJOR,N+1,1,M,N+1,ipiv,RHS,N+1);
    if (info != 0)
      throw std::runtime_error("fitCurve<2,4,ppForm>() - DGESV");
    //assemble the spline
    res.knots_ = t;
    res.polys_.resize(N);
    for (int i = 0 ; i < N ; i++){
      Vec<double,1> K({(v[i+1][1]-v[i][1])/(t[i+1]-t[i])});
      Vec<double,1> m1({RHS[i]});
      Vec<double,1> m2({RHS[i+1]});
      res.polys_[i][0] = (m1 + m2 - (K * 2))/((t[i+1] - t[i])*(t[i+1] -
                                   t[i]));
      res.polys_[i][1] = (K * 3 - (m1 * 2) - m2)/(t[i+1] - t[i]);
      res.polys_[i][2] = m1;
      res.polys_[i][3] = Vec<double,1> {v[i][1]};
    }
    //release memory
    delete [] M;
    delete [] RHS;
    delete [] ipiv;
  }
  return res;
}

template<>
Spline<1,2,ppForm> interpolate(const InterpConditions<double,1>& IC,
BCType type, std::pair<int,int> infopair){
  const std::vector<std::vector<double> >& Table = IC.getTable();
  std::vector<Vec<double,2> > Data;
  for (int i = 0 ; i < (int)Table[0].size() ; i++)
    Data.push_back(Vec<double,2>({Table[0][i],Table[1][i]}));
  return fitCurve<1,2>(Data);
}

template<>
Spline<1,4,ppForm> interpolate(const InterpConditions<double,1>& IC,
                               BCType type, std::pair<int,int> infopair){
  assert(type != complete && type != second);
  if (type == mixed)
    assert(infopair.first == 0 && infopair.second == 0);
  const std::vector<std::vector<double> >& Table = IC.getTable();
  std::vector<Vec<double,2> > Data;
  for (int i = 0 ; i < (int)Table[0].size() ; i++)
    Data.push_back(Vec<double,2>({Table[0][i],Table[1][i]}));
  return fitCurve<1,4>(Data,type);
}

template<>
Spline<1,4,ppForm> interpolate(const InterpConditions<double,2>& IC,
                               BCType type, std::pair<int,int> infopair){
  assert(type != second);
  if (type == mixed)
    assert(infopair.first != 2 && infopair.second != 2);
  const std::vector<std::vector<double> >& Table = IC.getTable();
  const std::vector<std::vector<bool> >& BoolTable = IC.getBoolTable();
  std::vector<Vec<double,2> > Data;
  const int Size = (int)Table[0].size();
  if (type == complete || (type == mixed && infopair.first == 1))
    assert(BoolTable[2][0] == 1);
  if (type == complete || (type == mixed && infopair.second == 1))
    assert(BoolTable[2][Size-1] == 1);
  for (int i = 0 ; i < Size ; i++)
    Data.push_back(Vec<double,2>({Table[0][i],Table[1][i]}));
  const double start = Table[2][0];
  const double end = Table[2][Size-1];
  return fitCurve<1,4>(Data,type,start,end,infopair);
}

template<>
Spline<1,4,ppForm> interpolate(const InterpConditions<double,3>& IC,
                               BCType type, std::pair<int,int> infopair){
  const std::vector<std::vector<double> >& Table = IC.getTable();
  const std::vector<std::vector<bool> >& BoolTable = IC.getBoolTable();
  std::vector<Vec<double,2> > Data;
  const int Size = (int)Table[0].size();
  double start = 0; 
  double end = 0;
  if (type == complete || (type == mixed && infopair.first == 1)){
    assert(BoolTable[2][0] == 1);
    start = Table[2][0];
  }
  if (type == complete || (type == mixed && infopair.second == 1)){
    assert(BoolTable[2][Size-1] == 1);
    end = Table[2][Size-1];
  }
  if (type == second || (type == mixed && infopair.first == 2)){
    assert(BoolTable[3][0] == 1);
    start = Table[3][0];
  }
    
  if (type == second || (type == mixed && infopair.second == 2)){
    assert(BoolTable[3][Size-1] == 1);
    end = Table[3][Size-1];
  }
    
  for (int i = 0 ; i < Size ; i++)
    Data.push_back(Vec<double,2>({Table[0][i],Table[1][i]}));
  
  
  return fitCurve<1,4>(Data,type,start,end,infopair);
}


//ppForm end
//------------------------------------
//BForm begin


template<int Order>
class Spline<1,Order,BForm>{
public:
  using Bvec = std::vector<Bspline<1,Order> >;
protected:
  Bvec basis_;
  std::vector<double> coes_;
public:
  Spline() = default;
  Spline(const Bvec& bs, const std::vector<double>& coes);
  ~Spline() = default;
  void setbasis(const Bvec& bs);
  void setcoes(const std::vector<double>& coes);
  const Bvec& getbasis() const;
  const std::vector<double>& getcoes() const;
  const double operator()(double t) const;
  void show(std::ostream& os = std::cout) const;
  void addbasis(const std::vector<double>& knts);
  template<int Ord>
  friend Spline<1,Ord,BForm> fitCurve(const std::vector<double>& knots,
                                      const std::vector<double>& x,
                                      const std::vector<double>& y);
  template<int Ord>
  friend Spline<1,Ord,BForm> fitCurve(const std::vector<double>& x,
                                      const std::vector<double>&
                                      y);
  template<int Ord>
  friend void fnplt(const Spline<1,Ord,BForm>& Bs, const std::string
                    &_file, const int num, const bool overwrite);
  
};

template<int Order>
inline Spline<1,Order,BForm>::Spline(const Bvec& bs, const std::vector<double>& coes):basis_(bs),coes_(coes){}

template<int Order>
inline void Spline<1,Order,BForm>::setbasis(const Bvec& bs){
  basis_ = bs;
}

template<int Order>
inline void Spline<1,Order,BForm>::setcoes(const std::vector<double>& coes){
  coes_ = coes;
}

template<int Order>
inline const std::vector<Bspline<1,Order> >& Spline<1,Order,BForm>::getbasis() const{
  return basis_;
}

template<int Order>
inline const std::vector<double>& Spline<1,Order,BForm>::getcoes() const{
  return coes_;
}

template<int Order>
const double Spline<1,Order,BForm>::operator()(double t) const{
  const int num = coes_.size();
  double sum = 0;
  for (int i = 0 ; i < num ; i++)
    sum += coes_[i]*(basis_[i])(t);
  return sum;
}

template<int Order>
void Spline<1,Order,BForm>::show(std::ostream& os) const{
  const int num = coes_.size();
  for (int i = 0 ; i < num ; i++){
    basis_[i].show(os);
    os << "with coefficient = " << coes_[i] << endl;
  }
}

template<int Order>
void Spline<1,Order,BForm>::addbasis(const std::vector<double>& knts){
  const int numknots = knts.size();
  assert(basis_.empty() && numknots > Order);
  const int size = numknots - Order;
  std::vector<double> knots;
  for (auto it = knts.cbegin();it != knts.cend(); it++){
    knots.push_back(*it);
  }
  std::sort(knots.begin(),knots.end()); 
  basis_.resize(size);
  for (int i = 0 ; i < size ; i++){
    for (int j = 0 ; j <= Order ; j++){
      basis_[i].insert(knots[i+j]);
    }
  }
  int k = numknots - 1;
  std::vector<Bspline<1,Order> > vB1(k);
  std::vector<bool> flag(k);
  for (int i = 0 ; i < k ; i++ ){
    if (!equal(knots[i],knots[i+1])){
      std::set<double> tmpset;
      tmpset.insert(knots[i]);
      tmpset.insert(knots[i+1]);
      vB1[i].setknots(tmpset);
      std::vector<Polynomial<Order,double> > tmppoly(1);
      Polynomial<Order,double> p1(1);
      tmppoly[0] = p1;
      vB1[i].setpolys(tmppoly);
    }
    else{
      std::set<double> tmpset;
      tmpset.insert(knots[i]);
      vB1[i].setknots(tmpset);
      flag[i] = true;
    }
  }
  std::vector<std::vector<Bspline<1,Order> > > table(k);
  table[k-1] = vB1;
  for (k-- ; k >= size ; k--){
    std::vector<Bspline<1,Order> > vB(k);
    for (int i = 0 ; i < k ; i++){
      if (table[k][i].empty() == 0 || table[k][i+1].empty() == 0){
        std::set<double> tmpset = table[k][i].getknotsCopy();
        tmpset.insert(knots[i+numknots-k]);
        vB[i].setknots(tmpset);
        std::vector<Polynomial<Order,double> > tmppoly(numknots-k);
        if (table[k][i].empty()){
          std::vector<Polynomial<Order,double> > vp2 =
          table[k][i+1].getpolysCopy();
          Polynomial<2,double>
      tmp2({-1.0/(knots[i+numknots-k]-knots[i+1]),knots[i+numknots-k]/(knots[i+numknots-k]-knots[i+1])});
          for (int j = 0 ; j < numknots - k - 1; j++)
            vp2[j] = adjust<Order,Order+1>(vp2[j]*tmp2);
          tmppoly[0] = Polynomial<Order,double>(0);
          for (int j = 1 ; j < numknots-k-1; j++)
            tmppoly[j] = vp2[j-1];
          tmppoly[numknots-k-1] = vp2[numknots-k-2];
          vB[i].setpolys(tmppoly);
        }
        else if (table[k][i+1].empty()){
          std::vector<Polynomial<Order,double> > vp1 =
          table[k][i].getpolysCopy();
          Polynomial<2,double>
        tmp1({1.0/(knots[i+numknots-k-1]-knots[i]),-knots[i]/(knots[i+numknots-k-1]-knots[i])});
          for (int j = 0 ; j < numknots - k - 1; j++)
            vp1[j] = adjust<Order,Order+1>(vp1[j]*tmp1);
          tmppoly[0] = vp1[0];
          for (int j = 1 ; j < numknots-k-1; j++)
            tmppoly[j] = vp1[j];
          tmppoly[numknots-k-1] = Polynomial<Order,double>(0);
          vB[i].setpolys(tmppoly);
        }
        else{
          std::vector<Polynomial<Order,double> > vp1 = table[k][i].getpolysCopy();
          std::vector<Polynomial<Order,double> > vp2 = table[k][i+1].getpolysCopy();
          Polynomial<2,double>
        tmp1({1.0/(knots[i+numknots-k-1]-knots[i]),-knots[i]/(knots[i+numknots-k-1]-knots[i])});
          Polynomial<2,double>
      tmp2({-1.0/(knots[i+numknots-k]-knots[i+1]),knots[i+numknots-k]/(knots[i+numknots-k]-knots[i+1])});
          for (int j = 0 ; j < numknots - k - 1; j++){
            vp1[j] = adjust<Order,Order+1>(vp1[j]*tmp1);
            vp2[j] = adjust<Order,Order+1>(vp2[j]*tmp2);
          }
          tmppoly[0] = vp1[0];
          for (int j = 1 ; j < numknots-k-1; j++)
            tmppoly[j] = vp1[j] + vp2[j-1];
          tmppoly[numknots-k-1] = vp2[numknots-k-2];
          vB[i].setpolys(tmppoly);
        }
      
      }
      
    }
    table[k-1] = vB;
  }
  for (int j = 0 ; j < size ; j++){
    const std::vector<Polynomial<Order,double> >& respoly =
    table[k][j].getpolys();
    basis_[j].setpolys(respoly);
    auto it = basis_[j].getpolys().cbegin();
    for (int i = j ; i < j + Order ; i++ ){
      if (flag[i]){
        basis_[j].erase(it);
      }
      else
        it++;
    }
  }
  
}

template<int Ord>
Spline<1,Ord,BForm> fitCurve(const std::vector<double>& knots,
                                      const std::vector<double>& x,
                                      const std::vector<double>& y){
  const int numknots = knots.size();
  const int xknots = x.size();
  const int yknots = y.size();
  assert(numknots - xknots == Ord && xknots == yknots);
  Spline<1,Ord,BForm> res;
  res.addbasis(knots);
  const std::vector<Bspline<1,Ord> >& basis = res.getbasis(); 
  //prepare matrix
  double* M = new double[xknots*xknots];
  int dflag = 0;
  for (int i = 0 ; i < xknots ; i++){
    for (int j = 0 ; j < xknots ; j++){
      M[i+j*xknots] = basis[j].rder(x[i],dflag);
    }
    if (i < xknots -1 && equal(x[i],x[i+1])){
      dflag++;
    }
    else 
      dflag = 0;
    if (dflag > 0){
      const int r = std::count(knots.begin(),knots.end(),x[i]);
      assert(Ord - 1 - r - dflag >= 0);
    } 
  }
  // prepare RHS
  double* RHS = new double[yknots];
  for (int i = 0 ; i < yknots ; i++)
    RHS[i] = y[i];
  //solve the equation
  int* ipiv = new int[xknots];
  auto info = LAPACKE_dgesv(LAPACK_COL_MAJOR,xknots,1,M,xknots,ipiv,RHS,yknots);
  if (info != 0)
    throw std::runtime_error("fitCurve<1,Ord,BForm>() - DGESV");
  std::vector<double> coe(RHS,RHS+yknots);
  res.setcoes(coe);
  delete [] M;
  delete [] RHS;
  delete [] ipiv;
  return res;
}

template<int Ord>
Spline<1,Ord,BForm> fitCurve(const std::vector<double>& x,
                             const std::vector<double>& y){
  std::vector<double> knots = aptknt(x,Ord);
  return fitCurve<Ord>(knots,x,y);
}

template<int Ord>
void fnplt(const Spline<1,Ord,BForm>& Bs, const std::string
           &_file, const int num = 200, const bool overwrite = 0){
  const std::vector<Bspline<1,Ord> >& basis = Bs.getbasis();
  const double begin = *(basis.front().getknots().cbegin());
  const double end = *(--(basis.back().getknots().cend()));
  const double h = (end - begin)/num;
  double t = begin;
  std::vector<double> y(num+1);
  std::ofstream os;
  if (overwrite)
    os.open(_file,ios::app);
  else
    os.open(_file);
  os << "t=[";
  for (int i = 0 ; i < num ; i++){
    y[i] = Bs(t);
    os << t << ",";
    t += h;
  }
  os << t;
  y[num] = Bs(end);
  os << "];\n";
  os << "y=[";
  for (int i = 0 ; i <= num ; i++){
    os << y[i] << ",";
  }
  os << "];\n" << "plot(t,y);";
  os.close();
}

//BForm end
//-------------------------------
//cardinalB begin

template<int Order>
class Spline<1,Order,cardinalB>: public Spline<1,Order,BForm>{
public:
  using BaseClass = Spline<1,Order,BForm>;
  using Bvec = std::vector<Bspline<1,Order> >;
public:
  Spline() = default;
  Spline(const Bvec& bs, const std::vector<double>& coes);
  ~Spline() = default;
  template<int Ord>
  friend Spline<1,Ord,cardinalB> fitCurve(const std::vector<int>& knots,
                                      const std::vector<double>& x,
                                      const std::vector<double>& y);
  template<int Ord>
  friend Spline<1,Ord,cardinalB> fitCurve(const int& start, const int& end,
                                      const std::vector<double>& x,
                                      const std::vector<double>& y);
};

template<int Order>
Spline<1,Order,cardinalB>::Spline(const Bvec& bs, const std::vector<double>& coes):BaseClass(bs,coes){};

template<int Ord>
Spline<1,Ord,cardinalB> fitCurve(const std::vector<int>& knots,
                                 const std::vector<double>& x,
                                 const std::vector<double>& y){
  std::vector<double> nknots(knots.begin(),knots.end());
  Spline<1,Ord,BForm> tmp = fitCurve<Ord>(nknots,x,y);
  Spline<1,Ord,cardinalB> res(tmp.getbasis(),tmp.getcoes());
  return res;
}

template<int Ord>
Spline<1,Ord,cardinalB> fitCurve(const int& start, const int& end,
                                 const std::vector<double>& x,
                                 const std::vector<double>& y){
  const int num = x.size();
  assert(end - start + 1 == num + Ord);
  std::vector<double> nknots;
  for (int i = start; i <= end ; i++)
    nknots.push_back(i);
  Spline<1,Ord,BForm> tmp = fitCurve<Ord>(nknots,x,y);
  Spline<1,Ord,cardinalB> res(tmp.getbasis(),tmp.getcoes());
  return res;
}

#endif // _SPLINE_H_
